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ABSTRACT 


Flame  initiation  and  propagation  through  an  air-fuel  vapor-fuel  drop 
system  is  numerically  modeled  in  a  cylindrical  one-dimensional  closed  com¬ 
bustor.  An  unsteady  formulation  of  the  flow  problem  eliminates  the  cold¬ 
boundary  difficulty  and  gas  phase  ignition  problem.  A  velocity  lag  between 
the  gas  and  the  liquid  phase  is  allowed  and  unsteady  heat  transfer  to  the 
droplets  is  taken  into  account.  The  surface  temperature  of  the  droplet  is 
evaluated  by  using  an  unsteady  spherically  symmetric  formulation  of  the 
droplet  heat  conduction  problem  with  no  internal  motion  and  with  a  time- 
varying  heat  flux  specified  at  the  surface  as  a  boundary  condition.  Results 
have  been  obtained  for  two  commercially  important  fuels,  namely,  n-octane 
and  n-decane.  The  activation  energy  and  the  preexponential  factor  in  the 
Arrhenius  type  expression  for  chemical  rate,  along  with  initial  temperature, 
initial  droplet  size,  stoichiometric  ratio,  and  diffusivity  are  parametrically 
varied  and  flame  speed  and  flame  temperatures  are  observed.  Flame  speed  is 
seen  to  increase  with  increasing  preexponential  factor,  with  decreasing  ac¬ 
tivation  energy,  with  increasing  ambient  temperature,  with  decreasing  initial 
droplet  radii,  and  with  increasing  diffusivities .  It  is  also  observed  that 
unlike  premixed  combustion,  heterogeneous  combustion  gives  rise  to  local 
variations  in  the  stoichiometric  fuel-to-air  ratio  in  the  axial  direction. 

This  phenomenon  could  give  rise  to  a  secondary  diffusion  flame  in  the  wake 
of  a  propagating  flame  and  produce  local  variations  in  the  flame  temperature. 
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I.  INTRODUCTION 


Modeling  effort  in  spray  combustion  has  received  attention  only  in 
the  past  two  or  three  decades.  Although  Wolfhard  and  Parker  [1]  reported  a 
qualitative  study  of  kerosene  spray  flame  propagation  in  1949,  the  first 
systematic  experimental  study  was  undertaken  by  Burgoyne  and  Cohen  [2]  in 
1954.  They  used  small  (up  to  10  microns  diameter)  tetralin  droplets  to 
study  the  effect  of  droplet  sice  on  flame  speed  and  concluded  that  flame 
speed  increases  with  decreasing  droplet  sice.  Similar  observations  were 
reported  by  Micutani  and  Ogasawara  [3]  with  much  larger  droplets  (diameter 
of  the  order  of  100u).  Later  Micutani  and  co-workers  studied  the  effect  of 
adding  relatively  nonvolatile  fuel  droplets  to  volatile  fuel-air  mixtures 
in  open  burner  flames  [4],  in  closed  spherical  vessels  [5],  and  in  spark 
ignition  engines  [6].  In  each  of  these  situations,  the  lean  flammability 
limit  was  extended  by  such  addition,  provided  the  droplet  sizes  were  large 
enough.  For  smaller  droplets,  the  behavior  shifted  towards  the  premixed 
behavior.  An  increase  in  flame  speed  was  also  reported  with  spray  addition. 
In  a  similar  study,  Polymeropoulos  and  Das  [7],  by  photographic  observation 
of  an  inverted  cone  flame  using  kerosene  air  spray,  concluded  that  the 
degree  of  atomization  of  spray  first  increases  the  flame  speed  and  then 
decreases  it  to  the  premixed  combustion  value. 

By  using  the  principle  of  the  Wilson  Cloud  Chamber,  Hayashi  and  Kumagai 
recently  obtained  a  uniformly  distributed  monodisperse  spray  of  ethanol  and 
n-octane.  In  their  earlier  study  with  this  technique  [8],  using  small 
ethanol  droplets,  they  observed  similar  effects  of  droplet  size  on  flame 
speed  as  observed  in  the  earlier  studies.  In  their  later  study  [9],  with 
droplets  of  up  to  40  microns  diameter,  they  observed  that  the  presence  of 
liquid  phase  in  lean  mixtures  resulted  in  flame  speeds  smaller  than  the 


2 

premixed  flame  speed  at  the  same  overall  fuel  mole  fraction.  The  rich 
mixtures,  on  the  other  hand,  showed  an  increase  in  flame  speed  compared 
to  the  premixed  flame  speed. 

Initial  analytic  modeling  was  limited  to  rather  simplified  cases 
due  to  lack  of  in-depth  understanding  of  the  related  processes.  Thus, 

Williams  [10]  assumed  a  steady  state,  monodisperse  spray  without  a  velocity 
lag  between  the  droplets  and  the  gas  and  analyzed  one-dimensional  spray 
deflagration.  He  further  simplified  his  analysis  by  assuming  the  fuel 
droplets  to  be  either  completely  vaporized  ahead  of  the  flame  or  neg¬ 
ligibly  vaporized.  In  the  first  case,  the  flame  propagation  was  essen¬ 
tially  premixed,  while  in  the  second  case  a  relay  mechanism,  in  which  the 
flame  propagated  from  one  droplet  to  another,  was  proposed.  Inverse  pro¬ 
portionality  between  droplet  size  and  flame  speed  was  obtained  for  a 
fixed  amount  of  fuel  in  the  second  case.  Mizutani  [11]  performed  a  much 
more  realistic  study  by  considering  a  polydisperse  spray,  turbulent  dif- 
fusivities,  and  velocity  lag  between  the  gas  and  the  droplets.  Droplet 
temperature  was  allowed  to  vary  temporally.  However,  prevaporization 
was  not  allowed  and  a  rigorous  criterion  for  droplet  ignition  was  not 
used;  ignition  occurred  whenever  the  droplet  temperature  exceeded  a 
critical  value.  Flame  speed  was  seen  to  increase  with  decrease  in  droplet 
size,  with  increase  in  turbulent  intensity,  and  with  increase  in  initial 
temperature.  Polymeropoulos  [12]  analyzed  one-dimensional  steady-state, 
laminar  flame  propagation  in  a  fuel  spray  system.  Prevaporization  as  well 
as  temporal  heating  of  droplets  was  allowed.  He  used  the  concept  of  igni¬ 
tion  temperature  for  the  droplet  ignition  and  again  showed  am  increase  in 
flame  speed  with  decrease  in  the  droplet  diameter  up  to  a  droplet  size 
beyond  which  the  flame  speed  was  seen  to  decrease  to  the  premixed  value. 


Bracco  [13]  studied  the  phenomenon  of  flame  initiation  and  ignition  by 
treating  the  gas  pha-e  as  unsteady.  He  modeled  fuel  droplets  falling 
under  gravity  through  a  stagnant  layer  of  hot  air.  The  droplets,  however, 
were  considered  to  be  at  constant  temperature  and  undergoing  pure  vaporiza¬ 
tion.  The  droplet  ignition  problem  was  argued  to  be  unimportant.  Computa¬ 
tions  were  conducted  for  cases  in  which  (a)  the  droplets  were  completely 
prevaporized  before  reaching  the  hot  air  zone,  (b)  were  partially  pre¬ 
vaporized,  and  (c)  were  ignited. 

In  this  report,  a  numerical  study  is  presented  for  the  combustion  of 
a  single-component,  polydisperse  fuel  spray  in  a  one-dimensional  closed 
combustor  with  laminar  flow.  The  spray  is  assumed  to  be  dilute  so  that 
interaction  between  the  droplets  is  negligible.  Transient  heating  and 
pre-vaporization  of  the  droplets  is  allowed.  The  coupling  between  droplet 
vaporization  and  environmental  conditions,  shown  to  be  significant  in  a 
cold  environment  (Law  [14]),  is  considered.  Lag  between  gas  and  droplet 
velocities  is  allowed  and  the  droplet  motion  is  described  by  a  drag  law. 

The  effect  of  convection  is  described  by  a  Ranz-Marshal 1  type  correlation. 
The  need  for  a  gas-phase  ignition  criterion  has  been  eliminated  by  an  un¬ 
steady  formulation  of  the  gas  phase.  The  droplet  ignition  criterion  has 
been  based  on  the  concept  of  ignition  Damkohler  number,  as  derived  by 
Law  [15].  Using  the  above  model,  the  phenomena  of  ignition,  flame  initia¬ 
tion  and  propagation  have  been  studied.  In  addition,  the  effects  of  initial 
temperature,  initial  droplet  size,  equivalence  ratio,  fuel  volatility,  and 
activation  energy,  on  these  phenomena  have  been  observed. 

II.  MATHEMATICAL  FORMULATION 

Firstly,  the  important  assumptions  used  in  the  present  model  will  be 


listed.  Semi -empirical  correlations  are  used  in  the  expressions  for  drop¬ 
let  drag  law  and  droplet  vaporization  rate.  It  is  assumed  that  the  vis¬ 
cous  dissipation  rate  is  negligible  and  the  kinetic  energy  of  the  flow  is 
negligibly  small  as  compared  to  its  thermal  energy.  As  a  result,  the 
pressure  is  assumed  to  be  spatially  constant.  The  wall  effects  are  ignored, 
i.e.,  a  free  slip  adiabatic  wall  is  considered.  The  gas  phase  is  assumed 
to  be  thermally  and  calorically  perfect  with  equal  (averaged)  specific 
heats  for  each  species.  Binary  diffusion  coefficients  for  each  pair  of 
species  are  taken  to  be  equal,  thermal  diffusion  is  neglected,  and  Fick's 
law  of  diffusion  is  employed  with  the  diffusion  coefficient  taken  as  3 
function  of  temperature  and  pressure.  Also,  Fourier's  law  of  heat  conduc¬ 
tion  is  used  and  thermal  conductivity  is  expressed  in  terms  of  the  diffu¬ 
sion  coefficient  by  assuming  the  Lewis-Semenov  number  to  be  unity.  Radi¬ 
ative  heat  transfer  is  assumed  negligible.  Latent  heat  of  vaporization 
and  heat  of  combustion  are  assumed  constant.  The  Clasius-Clapeyron  relation 
is  used  to  determine  the  fuel  vapor  concentration  at  the  droplet  surface. 

The  gas  phase  chemical  reaction  is  assumed  to  be  governed  by  a  single-step 
second-order  kinetics  of  the  type; 

Fuel  +  Oxidizer  -►  Products 

Further,  in  the  analysis  of  the  liquid  phase,  the  gas  film  surrounding  the 
droplet  is  considered  as  quasi-steady.  The  polydisperse  spray  has  been 
approximated  by  superposition  of  mono-disperse  sprays. 

Two  basic  approaches  have  been  used  in  the  literature  to  formulate 
the  spray  combustion  problem.  The  first  approach  is  to  consider  a  steady- 
state  flow  problem  and  thus  solve  a  set  of  ordinary  differential  equations. 
The  formulation  is  faced  with  the  cold  boundary  difficulty,  which  is  usually 
removed  by  suppressing  the  chemical  reactions  and  droplet  vaporization  below 
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certain  temperature.  The  second  approach,  adopted  in  the  present  formula¬ 
tion,  treats  the  whole  gas  phase  problem  as  unsteady.  This  approach  has 
wider  application  and  is  easier  to  handle  numerically. 


Gas-Phase  Flow  Equations 


Considering  the  source  terms  due  to  evaporation  and  heat  transfer  from 
the  droplets,  and  using  the  assumptions  listed  above,  we  can 
obtain  the  following  conservation  equations  (Reference  16)  for  the  overall 
mass  density  (p),  species  mass  fractions  (Y^ ,  momentum  and  energy  re¬ 
spectively: 


3p  .3  ,  ,  • 

It  +  (pu)  =  nmv  =  E  nimv 

i  =  l  i 

3(PYi)  3  a  3Y. 

3  +  (puyi)  "  (pd  T~)  =  5.  nm  +  W 

®  ax  i  3x  3x  iF  v  chem. 


where  i  corresponds  to  fuel,  oxidizer  and  product  respectively. 


|£=  o 

3x 


(1) 

(2) 

(3) 


pCn  If  +  pC  u  I1  '  f-  (pDC  — )  - 
P  31  p  3x  3x  p  3xJ  dt 

N 

=  W  Q  -  E  n.m  [C  (T-T  )  +  L!] 
f  i=l  1  l  p  s-  iJ 

Also  we  have  the  ideal  gas  equation 


5iF  in  Equation  (2)  is  the  kronecker  delta  and  is  1  when  i  equals  F  and 
is  0  otherwise.  I!  in  Equation  (4)  is  the  effective  latent  heat  of  vapor 
ization  and  includes  unsteady  conduction  of  he3t  to  the  droplets.  The 
source  terms  in  Equations  (1)  to  (4)  arc  given  as 


(6) 


1-Y  (l-5fW6fY 

(1  ♦  .  3/Re~ )  4-nR.  PD  in( - — - - — ) 

1  1  1'YFs. 

l 


L! 


(l-VFs  )(C  (T-T  ).VS  Y  (}) 

_ 1  r _ 1 _ 

Y_  +  v6.Y  -  Y_(l-6-) 
Fs^  f  o  F  f 


(7) 


WF  = 


A  pV/^ 


(8) 


o 

Both  m  and  L!  are  obtained  (Reference  17)  through  a  quasi-steady  analysis 
vi 

of  gas  phase  surrounding  a  droplet.  6^  =  0  corresponds  to  the  pure  vaporiza- 

tion  case,  whereas  6^-1  corresponds  to  the  droplet  burning  case.  The  value 

of  is  determined  by  a  droplet  ignition  criterion,  as  obtained  by  Law  [15]. 

Yps  ,  the  fuel  vapor  mass  fraction  at  the  droplet  surface  is  related  to  T 
si  si 

through  the  Clasius-Clapeyron  relation. 

An  additional  equation  for  pressure  is  obtained  by  integrating  the 
energy  equation.  Using  the  facts  that  pressure  is  a  function  of  time 
only,  velocities  vanish  at  the  closed  ends  of  the  combustor  and  that  there 
is  no  heat  loss  from  the  ends  of  the  combustor,  Equation  (4)  can  be  inte¬ 
grated  to  yield 


l  N 


S-^r-w/  V*-J  ,z,  VY(Li  ■  Ts.)dx) 

* o  *o  1=1  1  1 


(9) 


Similarly,  integration  of  the  overall  continuity  equation  yields  the  ex¬ 
pression  for  the  gas  velocity 

x  ,  rx 

nm  dx  -  -rr  f 
'0  *0 

It  was  anticipated  by  studying  the  nature  of  the  first  integral  on 


i  t  x  a  rx 
uCx)  =  ^00  {J  ™vdx  '  at  J  pdx} 


(10) 


the  right  hand  side  of  Equation  (9)  and  later  verified  by  computer  results 
that  considerable  error  could  be  caused  in  the  numerical  evaluation  of  this 
term  [18] . 
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A  better  approach  is  to  use  Equation  (10)  at  x  =  l,  i.e., 
i  l 

u(z)  =  fw  {/o  n\dx  =  If  /o  pdxl  =  0  (11) 

The  above  equation  has  been  equated  to  zero  since  the  velocity  is  identi¬ 
cally  equal  to  zero  at  the  closed  end  (x  =  1).  The  gas  density  p  is  de¬ 
termined  if  the  gas  temperature,  the  pressure  and  the  gas  composition  is 
known  through  the  perfect  gas  equation.  The  temperature  and  the  composition 
are  known  from  the  solutions  of  the  energy  and  species  continuity  equations, 
respectively.  Equation  (11)  can  thus  be  used  to  evaluate  pressure.  It  may 
be  noted  that  for  a  combustor  open  at  x  =  l,  the  pressure  becomes  a  constant 
and  then  u(t)  can  be  obtained  by  using  Equation  (10). 

Using  the  transformation: 


1-v  r 
K  =  Tp  ——  =  Tp 


(12) 


the  energy  equation  (4)  is  rewritten  as 
3K 


3K  +  _ 

3t  U  3x 


W  Q  N  nimv 


(13) 


This  eliminates  the  pressure  term  in  the  energy  equation. 


Liquid-Phase  Equations 

For  each  group  of  monodisperse  droplets,  the  general  spray  distribu¬ 
tion  function,  which  may  be  a  function  of  time,  distance,  velocity,  and 
radii  of  the  droplets,  can  be  separated  into  three  simpler  equations.  These 
three  equations,  expressing  the  conservation  of  droplet  number  density, 
droplet  mass  and  droplet  momentum,  are  as  follows: 


3n.  3(n.u.) 

+  =  o, 

3t  3x 


3R.  3R.  v. 

— i  +  u  — i  -  x- 

3t  l.  3x 
l 


4lrRip* 


i  =  1,2,.  ...N 

i  =  1,2,. ...N 


(14) 

(15) 
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3u  . 
_ ii 

3t 


3u 


+  u 


Zi 


i  i  3x 


=  Fi, 


i  =  1,2,. ...N 


(16) 


where 


_.  81 

Fl  =  T 


(Rei)'Vi 


(u-ua)|u-u2i| 


(17) 


Re^  is  the  Reynolds  number  based  on  the  gas  kinematic  viscosity,  the  droplet 

relative  velocity  and  the  droplet  size. 

Assuming  that  T  is  known  (discussed  in  the  following  subsection), 
i 

we  then  have  (3N+S)  variables,  i.e.,  3N  spray  variables,  n^,  R^  and  U  ^ , 
three  species  mass  fractions  Yp,  Yq  and  Y  ,  gas  velocity  u,  gas  temperature  T, 
gas  pressure  p,  gas  density  p,  and  thermodynamic  variable  K.  Equations  (2),  (5), 
(10)  through  (16)  give  a  system  of  (3N+8)  equations  for  these  unknowns. 

The  initial  profiles  of  all  variables  are  taken  to  be  uniform  in  the  com¬ 
bustor.  Further,  the  initial  gas  and  droplet  velocities  are  assumed  to 
be  zero.  The  boundary  conditions  are  provided  as  follows: 

At  x  =  0 

3Y, 


3x 


•w 

o  =  _n  =  ii=u=u  =0 

3x  3x  3x  U  2.i 


n.  =  n. 

1  10 


*  m 
.  t  v. 


m  »-  v  . 

R.  =  R.  -  f  - y 

1  10  *0  4ttR~p 

u  Jl 


dt 


(18) 


At  x 


3x 


!!o  =  !!n  =  3K  =  0 

3x  3x  3x 


Droplet  Surface  Temperature 


Accurate  determination  of  the  droplet  surface  temperature  under  vapor¬ 
izing  and  combusting  conditions  is  very  important  in  spray  combustion  modeling. 


In  many  practical  applications,  droplets  undergo  a  period  of  pure  vapori¬ 
zation  before  ignition  takes  place.  The  instant  at  which  a  flame  is 
established  around  the  droplets  depends  on  m3ny  ambient  properties,  in¬ 
cluding  temperature  and  fuel  concentration.  The  amount  of  vapor  fuel 
depends  on  the  rate  at  which  a  droplet  vaporizes,  which  in  turn,  is  a 
strong  function  of  droplet  surface  temperature.  Since  many  of  the  droplet 
characteristics  change  significantly  upon  ignition,  it  is  obvious  that  the 
surface  temperature  histories  have  to  be  evaluated  accurately.  This 
makes  the  treatment  of  the  droplet  heat  transfer  problem  necessarily 
unsteady. 

Further,  in  considering  a  realistic  multicomponent  fuel  droplet, 
where  vaporization  rates  of  different  components  depend  on  their  distribu¬ 
tion  within  the  droplet,  determination  of  the  temperature  profiles  in  the 
interior  of  the  droplet  is  essential  [19,20,21].  Although  in  the  present 
study  we  are  considering  only  monocomponent  fuel  droplets,  the  extension 
of  the  present  model  to  the  case  of  multicomponent  fuel  droplets  is  of 
interest.  We  will  analyze  droplet  transient  heating  in  the  limit  of  pure 
conduction  where  variation  in  temperature  is  possible  in  the  interior  of 
the  droplet.  Results  for  the  more  realistic  case  of  axisynmetric  internal 
circulation  inside  the  droplet  have  been  obtained  recently  (Prakash  and 
Sirignano  [22,23]).  However,  these  results  were  not  available  at  the 
time  this  model  was  formulated  and  we  have  employed  the  results  of  Law 
and  Sirignano  [24]  and  Law  [25].  The  details  of  the  droplet  energy 
equation  with  appropriate  boundary  conditions  are  given  in  Reference  [18]. 
The  solution  of  that  equation  with  a  constant  heat  flux  at  the  surface  can 
be  found  in  Carslaw  and  Jaeger  [26].  For  time  varying  surface  heat  flux 
F(t)  an  analytical  solution  is  obtained  by  using  Duhamel's  Principle. 
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Thus,  for  the  surface  temperature  of  the  droplet,  the  following  expression 
is  obtained. 


...  ,;2  ,  7  -  sin  ra  „  .  p\ 

T(r,t)  =  -  4  2  “2 - -]F(t:D  ♦  3 f  F(x)dx 

r  n=l  «  sin  a 
n  n 


_  ®  sin  ra  rt  .  -a  (t-x) 

♦  l  1  - 2-  f  F(x)e  n  dx 

r  n=l  Sin  an  'o 


(19) 


where  an>  n  =  1,2,...,®  are  positive  roots  of  the  equation 


tan  a  =  a.  (20) 

2 

For  a  large  enough  number  Nq,  for  which  >>  1,  the  above  equation  can 

—  ‘  o 

be  written  with  finite  sum  as 

-  N  2  ' 

Xt  -  o  -t  „  ,  <*(t-x) 

F(x)dx  +2  Z  J  F(t)e  n  dx  (21) 

i  n=l'o 

where 

N 

0  2 

Cj  =  0.2  -  2  £  1/a^  (22) 

n=  1 

The  integrals  present  in  Equation  (21)  cannot  be  directly  evaluated 
numerically  in  convective  situations,  since  the  equation  has  been  derived 
in  coordinates  fixed  to  a  droplet  whereas  the  flow  equations  are  obtained 
in  the  coordinates  of  the  combustor  in  which  the  droplets  are  moving.  In 
a  numerical  scheme,  the  integrands  can  be  calculated  only  at  discrete  grid 
points  and  cannot  follow  a  droplet  unless  an  extensive  interpolation  scheme 
is  used  to  approximate  the  values  between  the  grid  points.  A  more  accurate 
approach  is  to  calculate  the  integrals  by  solving  the  differential  equations 
obtained  from  a  transformation  which  relates  the  two  coordinate  systems. 

We  define  new  variables  h(t)  and  g(t)  as 

-  •  r*  - 

h(t)  =  /  F(x)dx  (23) 

* o 


and 
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N  N  :  2 ,  . 

0  .  0  r  -  -an(t-'r) 

g(t)  i  Z  g^t)^  z  I  F(T)e  n  dr  (24) 

n=  1  n=  1  “*0 

Under  the  transformation  from  the  frame  of  reference  of  the  moving  droplet 

to  the  frame  of  reference  of  the  combustor,  the  variables  Ts(t),  F(t), 
h(t)  and  g(t)  transform  into  Ts(t,x),  F(T,x),  h(t,x)  and  g(t,x),  respec¬ 
tively.  The  following  differential  equations  for  h  and  g  can  be  immediately 
written 


3h  3h 

—  ♦  u  — 
3t  l  3x 


ai 

-4  F(t,x) 

R“ 


(25) 


!!ll  +  Ifn 
3t  Ui  3x 


l  2 

,2  an^n 


(26) 


Thus  the  integrals  in  Equations  (21)  can  be  evaluated  by  solving  equations 
(25)  and  (26).  The  boundary  conditions  for  these  equations  are  easily 
obtained,  realizing  that  u^ = 0  at  the  x  =  0  boundary,  by  direct  integra¬ 
tion  of  Equations  (25)  and  (26)  which  gives  h(0,t)  and  g^(0,t), 

n  =  l,2f...,N  .  Tg  is  then  given  by 

- - -  N 

-  -  o  2  ' 

Tg(t)  =  (0.2  -  2  Z  l/a^)F  +  3h  +  2g  (27) 

n=l 


III.  NWERICAL  CON'S  I  DEFLATIONS 

The  final  set  of  governing  equations  are  first  nondimensionalized  with 
respect  to  appropriate  characteristic  values  so  that  all  variables  are  of 
the  same  order  of  magnitude.  The  characteristic  values  are  based  on  the 
phenomenon  of  interest.  Thus,  the  characteristic  value  of  distance  is  based 
on  the  flame  thickness  and  the  characteristic  value  of  time  is  based  on  the 
values  of  diffusion,  vaporization,  and  chemical  times.  Using  the  variables 
in  the  ignition  zone,  these  are  estimated  to  be 


1 


'0(0.2  cm.) 


6 


T 


diffusion 


T 


T 


vaporization 
chemical  ^ 


~0(0.1  sec.) 
“0(0.01  sec.) 
-0(0.001  sec.) 


The  characteristic  values  of  other  variables  such  as  pressure,  temperature, 
density,  droplet  size,  etc. ,  are  conveniently  chosen  at  their  initial  ambient 
values  and  the  characteristic  values  of  derived  quantities  are  obtained  from 
dimensional  considerations. 

Since  the  governing  equations  constitute  a  set  of  nonlinear  partial 
differential  equations,  it  would  be  easier  to  solve  them  if  they  are 
linearized  in  some  way.  We  will  be  using  the  technique  of  quasi-linearization, 
first  introduced  by  Bellman  and  Kalaba  [27]  and  later  used  by  Lee  [28]  and 
Sharma  and  Sirignano  [29]. 

In  this  technique,  the  nonlinear  terms  are  expanded  in  a  Taylor  series 
with  respect  to  the  variable  and  its  derivatives  around  the  previously  known 
value.  Thus,  for  a  typical  non-linear  differential  equation 

3Y 

=  Z(Y,Y  ,Y  )  (28) 

3t  x  xx  v 

the  quasi- linearized  form  is 


(!r)k+1  =  zk  +  c|f)k(Yk+1  -  vk) 


♦  (H^)k(Yx+1  -  Yx3 


♦  (4p— )k(Yk+1  -  Yk  ) 
3Y  xx  xx 

XX 


(29) 


where  subscript  x  represents  the  partial  derivative  with  respect  to  x, 
superscript  k  refers  to  a  previously  known  value  and  (k+1)  refers  to  an 
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unknown  value.  This  linearization  process  has  an  added  advantage  of  gen¬ 
erating  a  set  of  uncoupled  differential  equations.  These  equations,  when 
approximated  by  an  implicit  finite  difference  method,  leads  to  a  tridiagonal 
matrix  as  opposed  to  the  linearization  technique,  used  by  Briley  and 
McDonald  [30],  which  leads  to  a  block  tridiagonal  system.  However,  the 
present  method  needs  iteration  to  obtain  the  solution  at  each  time  step. 

The  quasi-linearized  equations  were  replaced  by  finite-difference  ap¬ 
proximations  by  employing  an  implicit  scheme,  similar  to  Crank-Nicolson 
scheme.  However,  it  was  observed  that  the  scheme,  when  used  for  the  first 
spatial  derivative,  lead  to  an  oscillatory  behavior  for  the  variables 
governed  by  first  order  differential  equations.  A  satisfactory  analog  for 
the  first  derivative  was  found  to  be  as  follows: 


tsh.,.1/2  ■  ,Yl,i-r  Yl,i  *  •<“•“> 

where  the  subscripts  i  and  j  represent,  respectively,  the  iAx  and  jAt 
locations  in  space  and  time.  Then  the  finite-difference  approximation 
to  Equation  (29)  is 

CkYk+;  .  ,  ♦  AkYk+1  ,  +  BkYk+!  .  =  Dk  (31) 

l-l, j+1  i,j+l  i+l, J  +  l 


k+ 1 

where  Y  is  the  unknown  variable  which  is  being  solved  at  each  of  the 
k  k  k  k 

grid  points  and  A  ,  B  ,  C  and  D  are  coefficients  calculated  from  the 
known  values  cf  the  variables  at  the  previous  iteration.  If  ID  is  the 
total  number  of  grid  points  then  we  can  obtain  ID-2  such  equations  for 
i  =  2 , 3, . . . , ID- 1 .  Boundary  conditions  provide  two  more  equations.  Thus 
we  have  ID  equations  for  ID  unknowns  Y 

Flame  initiation  has  been  modeled  by  providing  a  zone  of  heated  gases 
at  one  end  of  the  combustor.  A  temperature  of  1500°k  and  a  length  of 


one  half  cm  for  the  hot  cone  has  been  found  sufficient  for  ignition  in 
most  cases.  A  corresponding  amount  of  fuel  has  been  removed  from  the  hot 
zone  so  as  to  maintain  the  adiabatic  flame  temperature  the  same.  To 
avoid  large  gradients  near  the  ignition  source,  the  initial  values  of  all 
variables  are  varied  linearly  from  the  ignition  source  values  to  ambient 
values  in  the  next  one-half  centimeter.  Then  using  Equation  (31)  with 
boundary  conditions  (18),  the  variables  Yp,  Y^,  Y^,  k,  n^,  and  u^ 

1 

are  sequentially  solved.  New  pressure  and  velocity  values  are  obtained 
from  Equations  (11)  and  (10),  and  T  is  obtained  from  Equations  (25), 

l 

(26)  and  (27).  Again,  an  implicit  numerical  scheme  is  used  to  solve 

s  t 

Equations  (25)  and  (26).  Once  all  variables  are  known  at  k+1  iteration, 
new  coefficients  for  Equation  (31)  are  obtained.  This  process  is  repeated 
until  convergence  is  achieved  at  each  of  the  grid  points  for  a  given  time. 
Then  time  is  incremented  and  the  same  process  is  repeated.  To  check  the 
accuracy  of  the  method,  the  value  of  Y^,  as  obtained  by  solving  equation 
for  Y^,  is  compared  with  the  value  of  the  same  as  obtained  from  relation 


The  grid  spacings  in  the  space  and  time  axes,  i.e.,  Ax  and  At,  have 
been  based  on  the  characteristic  values  of  length  and  time.  The  value  of 
At  varies  from  one-half  millisecond  to  10  microseconds  depending  on  the 
characteristic  value  of  the  chemical  and  evaporation  times  for  various  cases. 
Another  consideration  for  choosing  At  is  the  cost.  A  smaller  At  requires 
fewer  quasilinearization  iterations  (k-iterations)  at  any  time  and  vice 
versa.  Although  a  rigorous  criterion  for  determining  At  has  not  been  derived, 
experience  shows  that  modeling  is  most  economical  when  the  number  of  k- 
iterations  at  any  time  lies  between  2  to  5.  It  may  be  noted  here  that  in 
some  of  the  cases  when  time  steps  of  the  order  of  10  microseconds  must  be 


IS 


chosen,  complete  flame  propagation  in  the  combustor  may  require  as  many 

4 

as  10,000  time  iterations  (approximately  3  *  10  k- iteration) .  Therefore, 
in  many  cases  modeling  has  been  terminated  when  the  flame  has  propagated  a 
significant  portion  of  the  combustor  length. 

The  values  of  physio-chemical  constants  for  the  fuel  and  the  ambient 
gas  have  been  compiled  from  various  sources  [30,31,32].  Some  variation 
was  seen  in  the  values  of  the  latent  heat  of  vaporization  and  heat  of  com¬ 
bustion  of  the  fuel.  Therefore,  whenever  possible,  all  data  for  a  fuel 
has  been  collected  from  a  single  source.  The  averaged  constant  value  of 


0^  is  selected  to  yield  the  adiabatic  flame  temperature  upon  complete 
combustion  of  stoichiometric  fuel  mixture.  In  reality,  value  of  de¬ 
pends  on  the  gas  temperature  and,  to  a  lesser  degree,  on  the  gas  pressure 


and  composition. 


Therefore  a  constant  C 

P 


would  not  yield  accurate  results 


for  incomplete  combustion,  different  stoichiometric  ratios  or  different 


initial  conditions. 


Cases  with  various  values  of  the  preexponential  factor  A,  activation 
energy  E,  initial  temperature,  T^,  overall  equivalence  ratio  $,  diffusivity 
D,  droplet  size  R,  and  cases  of  homogeneous  combustion  have  been  modeled. 
Parametric  variation  of  A  and  E  is  important  since  the  value  of  these  para¬ 
meters,  for  the  fuels  under  consideration,  is  not  known  with  much 
confidence . 


IV.  RESULTS  AND  DISCUSSION 

Before  we  discuss  the  results,  it  may  be  indicated  that  after  all  the 
above  mentioned  cases  had  been  calculated,  some  en-ors  were  found  in  the 
computer  code.  In  addition,  few  modifications  were  incorporated  in  the  code. 


For  example,  the  finite  difference  scheme  for  the  liquid  phase  equations 
and  the  droplet  surface  temperature  equations  was  made  second  order  accurate 
This  resulted  in  a  reduction  of  observed  flame  speed  as  the  amount  of  numer¬ 
ical  diffusion  was  reduced.  Other  modifications  were  aimed  at  making  the 
code  more  efficient.  It  is  expected,  though  not  claimed,  that  all  these 
alterations  will  effect  only  the  quantative  nature  of  the  results  obtained 
from  the  old  code.  Consequently,  only  one  control  case  for  n-Decane  fuel 
droplets  of  40U  radius  has  been  repeated  with  the  new  code.  All  the  other 
results  and  conclusions  have  been  derived  from  the  old  code.  The  results 
of  the  repeated  case  are  shown  in  Figures  la  to  lh, 
where  the  profiles  of  fuel  main  functions,  oxidizer  mass  friction,  gas 
temperature,  gas  velocity,  droplet  velocity,  droplet  number  density, 
droplet  radius  and  the  reaction  rate  are  plotted  in  the  combustor  for 
several  time  values.  The  grid  size  for  this  control  run  was  .625  mm 
and  the  time  step  varied  from  10us  to  lOOys.  It  was  observed,  through 
trial  and  error,  that  the  characteristic  reaction  time  controlled  the  time 
step. 

Figures  la,  lb  and  lc  indicate  some  distinguishing  features  of  flame 
propagation  in  air- fuel  spray  mixtures.  Firstly,  as  seen  in  Figure  la  and 
lb,  the  mixture  equivalence  ratio,  which,  initially  is  less  than  unity 
everywhere  in  the  combustor,  can  greatly  vary  along  the  length  of  the  com¬ 
bustor  at  later  times:  Secondly,  the  present  results  indicate  the  presence 
of  a  secondary  diffusion  flame  in  the  wake  of  the  primary  propagating  flame. 
The  location  of  the  primary  flame,  a  pre-mixed  type  flame,  is  indicated  by  a 
sharp  peak  in  the  fuel  mass  fraction  profiles.  This  sharp  peak  is  caused  by 
the  suddenly  enhanced  vaporization  rate  in  the  vicinity  of  the  primary  flame 
The  location  of  the  secondary  flame,  a  diffusion  type  flame,  is  the  point 


where  both  the  fuel  and  oxidizer  mass  fractions  vanish.  This  location  is 
particularly  evident  in  the  profiles  for  t  >_  104  ms  in  Figures  la  and  lb. 
Also,  the  first  peak  in  the  temperature  profiles  for  t  >_  104  ms  indicates 
this  diffusion  flame.  The  above  features  can  be  attributed  to  the  rela¬ 
tive  motion  of  the  gas  and  the  droplets.  This  relative  motion  is  demon¬ 
strated  in  Figures  Id  and  le.  In  these  figures,  the  point  where  the  velocity 
reverses  sign  represents  the  location  of  the  primary  flame.  Thus,  the  flows 
of  gas  and  droplets  are  always  away  from  the  primary  flame.  However,  in 
case  of  the  secondary  diffusion  flame,  fuel  and  oxidizer  flow  towards  the 
flame. 

In  the  colder  or  upstream  section  of  the  combustor,  gas  temperature 
rises,  as  indicated  by  Figure  lc,  due  to  adiabatic  compression  of  the  gas. 
This  enhances  the  heat  transfer  rate  to  the  droplet.  As  a  consequence, 
the  droplet  surface  temperature  rises,  which,  in  turn,  increases  the  droplet 
vaporization  rate  in  the  colder  region.  In  addition,  because  the  droplets 
move  away  from  the  flame,  the  droplet  number  density  increases,  in  the 
colder  region  of  the  combustor  (see  Figure  If).  This  further  enhances 
the  vaporization  rate  in  that  region.  Consequently,  as  seen  in  Figure  lg, 
the  droplet  radius  decreases  at  a  continuously  increasing  Tate  in  that 
region . 

It  is  also  worth  mentioning  that  Figure  If  exhibits  an  oscillatory 
behavior  for  the  droplet  number  density  near  the  combustor  ends.  This  is 
acceptable  near  the  left  end  as  the  droplet  radius,  there,  is  negligibly 
small.  The  oscillatory  behavior  near  the  right  end  is  possibly  due  to 
some  inconsistency  in  using  the  boundary  conditions  for  the  liquid  phase 
equations.  This  is,  presently,  under  investigation.  Profiles  of  the 
surface  temperature  of  the  droplets  were  qualitatively  similar  to  the  gas 
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temperature  profiles.  The  maximum  value  of  droplet  temperature,  however, 
was  limited  by  the  boiling  point  of  the  fuel. 

The  reaction  rate  profiles,  shown  in  Figure  lh,  indicate  a  very 
narrow  reaction  tone  for  the  primary  flame.  Consequently,  the  time  step 
and  the  grid  sice  for  the  finite  difference  equations  were  controlled  by 
the  thickness  of  this  reaction  zone.  The  presence  of  the  secondary  flame, 
a  diffusion  type  flame  moving  in  a  direction  opposite  to  that  of  the  primary 
flame,  is  also  indicated,  in  Figure  lh,  by  the  non-zero  reaction  rates  near 
x  =  2  cm  and  for  t  =  160  and  580  ms.  The  reaction  rate  profiles  for  the 
primary  flame  also  demonstrate  the  unsteady  nature  of  flame  propagation  in 
sprays.  The  maximum  value  of  the  reaction  rate  shows  a  wide  variation  as 
the  primary  flame  propagates  through  the  combustor.  This  variation  is  re¬ 
lated  to  the  variation  in  the  local  equivalence  ratio  of  the  mixture  ($) 
in  the  vicinity  of  the  primary  flame.  Table  I  lists  these  $  values  along 
with  the  values  of  the  maximum  reaction  rate.  As  this  table  indicates, 
the  value  of  the  maximum  reaction  rate  rises  as  <p  approaches  unity.  This 
variation  in  the  reaction  rate  is  also  responsible  for  the  peculiar  nature 
of  the  temperature  profiles  between  the  two  temperature  peaks  corresponding 
to  the  primary  and  secondary  flames  (see  Figure  lc).  For  example,  for  the 
profile  at  t  =  580  ms,  the  low  temperature  point  between  the  two  peaks  is 
due  to  the  fact  that  the  maximum  value  of  reaction  rate  at  earlier  times, 
say  at  t  =  160  and  320  ms,  was  relatively  small.  The  unsteady  nature  of 
the  flame  propagation  is  also  manifested  in  Table  II,  which  gives  the 
flame  location,  the  fl3me  velocity  relative  to  the  combustor,  the  gas 
velocity  and  the  flame  front  velocity  relative  to  the  gas  for  several  t 
values.  All  these  values  pertain  to  the  primary  flame.  The  flame  location 
is  defined  as  the  point  where  the  gas  temperature  is  1000°k.  The  variation 


in  the  flame  velocity  relation  to  gas  is  related,  through  the  reaction 
rate,  to  the  variation  in  the  local  equivalence  ratio. 

In  the  rest  of  this  section,  the  effects  of  varying  A,  E,  Tq  and  D 
are  discussed.  The  results  of  two  other  runs,  one  with  a  pre-mixed  homo¬ 
geneous  mixture  and  the  other  with  two  mono-disperse  droplet  gToups,  are 
also  mentioned.  Since  all  of  these  results  were  obtained  by  using  the 
older  version  of  the  computer  code,  only  a  qualitative  discussion  can  be 
given  with  confidence. 

A  reduction  in  preexponential  factor  A  caused  a  reduction  in  the  ob¬ 
served  flame  velocities.  This  could  lead  to  an  extinguishment  of  the 
flame.  A  reduction  in  activation  energy  E  had  just  the  opposite  effect. 

In  addition,  the  secondary  diffusion  flame  would  appear  earlier  if  E  is 
reduced.  The  effect  of  increasing  the  initial  gas  temperature  T  ,  from 
300°k  to  600°,  is  shown,  3gain  qualitatively,  in  Figures  2a  and  2b.  These 
figures  indicate  that  the  variation  in  the  equivalence  ratio  is  not  quite 
as  severe  as  in  the  control  run.  The  reason  is  the  enhanced  heat  transfer 
rate  to  the  droplets,  which  increases  the  vaporization  rate  in  the  colder 
part  of  the  combustor.  In  other  words,  there  is  a  tendency  to  approach 
the  case  of  pre-mixed  homogeneous  type  of  combustion  away  from  the  igni¬ 
tion  source.  For  moderate  increase  in  T,  however,  we  still  have  a  fuel 
rich  region  near  the  ignition  source  and  an  oxidizer- rich  region,  formed 
behind  the  primary  flame  as  this  flame  travels  into  the  fuel-lean  region. 

A  diffusion  type  flame  appears  at  the  interface  of  these  two  regions.  The 
effects  of  increasing  mass  and  thermal  diffusion  coefficients  were,  as 
expected,  an  increased  flame  speed  and  an  increased  flame  thickness. 

Figure  3  gives  the  temperature  profiles  for  the  flame  propagation  in 
a  pre-mixed  homogeneous  mixture.  Again  the  results  are  only  qualitatively 
correct.  Important  features  are:  (11  As  expected,  there  is  no  secondary 


20 


diffusion  flame  and  (2)  the  primary  flame  tends  to  propagate  with  a  relatively 
steady  velocity. 

A  case  of  two  mono-disperse  droplet  groups,  which  represents 
the  simplest  poly-disperse  case,  was  also  modeled.  Initial  droplet 
sizes  of  4Cy  and  20u  were  used.  A  reduction  in  droplet  size  resulted 
in  an  increase  in  flame  speed.  Droplet  radius  and  velocity  profiles  are 
shown  in  Figures  4a  and  4b,  respectively.  Figure  4b  indicates  that  the 
smaller  droplets  can  follow  the  gas  more  easily. 


CONCLUSION'S 

The  phenomenon  of  one-dimensional  flame  propagation  in  air  and  fuel 
spray  mixtures,  with  transient  heating  of  the  droplets  allowed,  have  been 
successfully  modeled.  Some  preliminary  results  have  been  reported.  These 
results  are  more  useful  for  studying  the  processes  involved  in  a  qualita¬ 
tive  way  rather  than  for  making  quantitative  predictions.  It  is  clear 
that  the  process  is  essentially  unsteady  for  all  the  cases  considered 
with  the  exception  of  pre-mixed  case.  The  unsteadiness  is  due  to  (a)  the 
time-varying  stratification  caused  by  the  motion  of  the  droplets  and  (b)  the 
transient  heating  of  the  droplets. 

Before  obtaining  more  detailed  and  quantitatively  useful  results, 
further  refinement  and  optimization  of  the  computer  code  is  worth  consid¬ 
ering.  For  example,  a  non-iterative,  fully-implicit  scheme,  as  use  by 
Briley  and  McDonald  [30],  as  well  as  a  variable  grid  size  may  be  worth 
investigation.  In  future  work,  one  could: 

I.  Relax  the  non-essential  assumptions  like  constant  C^. 

II.  Include  the  effect  of  turbulence,  say,  by  using  a  k-e  type 


model . 


III.  Consider  certain  realistic  effects  such  as  (a)  addition  of 
gravity,  (b)  addition  of  wall  heat  transfer,  (c)  spark  ig¬ 
nition  process  and  (d)  a  combustor  open  at  both  ends,  useful 
for  studying  the  properties  of  steady  flames.  For  case  111(d), 
the  effect  of  droplet  diameter,  droplet  number  density,  over¬ 
all  equivalence  ratio  and  fuel  volatility  on  mechanism  of 
flame  propagation,  flame  speed  and  flammability  limits  can 
be  studied. 
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10  ms 


t  =  30  ms 


t  =  160  ms  t  =  520  ms  t  =  460  ms  t  *  580  ms  t  =  680  ms 


=  .031  W 


a  .08  W  =  .004  W  =  .004  W  =  .02  I  W 


1197 

3.9 
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.90 

480 
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870  .45 

642  .48 

519  .5 

435  .49 


1014  .40 

859  .46 

726  .50 

615  .55 

525  .58 


687  .72 

540  .77 

410  .75 


.053  W 


980  .93 

605  .93 

448  .85 

390  .8 


W  =  Maximum  value  of  the  Non-Dimensional  Reaction  Rate 
max 


Table  II.  Flame  Velocity  for  the  Control  Run 


985  .35 

579  .9 

441  .87 

396  .81 


x  =  Location  Flame  Velocity  Gas  Velocity  Flame  Velocity 
(t  in  ms)  f  T  i  nnn°t-  cm/sec  cm/sec  cm/sec 
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^cra  Combustor  Gas 
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